A Generalized Lambert Function of Two Arguments

The motivation for introducing this function comes from delay differential equations. The first order equation

$f^{'} (t) = a f (t - 1)$

with a unit displacement can be solved using the ordinary Lambert W function, defined implicitly as the inverse of the nonlinear transcendental equation

$W (x) e^{W (x)} = x$

Assuming an exponential solution of the form $e^{c t}$ , the differential equation simplifies to

$c = a e^{- c} \to c = W (a)$

Likewise, the first order equation

$f^{'} (t) = b f (t + 1)$

with a unit displacement in the opposite direction has a similar solution:

$c = b e^{c} \to c = - W (- b)$

What about an equation with displacements in both directions? With the same exponential solution, the equation

$f^{'} (t) = a f (t - 1) + b f (t + 1)$

simplifies to

$c = a e^{- c} + b e^{c}$

This equation is not soluble in terms of the ordinary Lambert function. Supposing there were some argument $c = W (x)$ that gave a solution to this equation, with the defining relation it would become

$W^{2} (x) (1 - \frac{a}{x}) = b x$

For $b = 0$ one can choose $x = a$ as expected, but for $a = 0$ one has

$W^{2} (- b) = - b^{2}$

which only holds for a finite number of points, not identically.

Since evaluation of the ordinary Lambert function occurs by numerical inversion in practice, one can simply define a new function as the solution to the equation

$W (x, y) = x e^{- W (x, y)} + y e^{W (x, y)}$

and treat it in the same fashion, taking into account the two special cases

$W (x, 0) = W (x) W (0, y) = - W (- y)$

There is also the interchange identity

$W (y, x) = - W (- x, - y)$

which is consistent with the two special cases.

One could, acting like a mathematician, change the sign of y in the defining equation to slightly simplify results here and later. Leaving the signs as given, however, helps to motivate expressions that are a bit more foolproof for purposes of explicit calculation. When it comes to accurate numerical evaluation it is better to act like a coder than a mathematician.

To have some idea of the nature of the problem in question, rewrite the defining equation as

$x e^{- W (x, y)} + y e^{W (x, y)} - W (x, y) = 0$

and visualize the real (blue) and imaginary (red) parts of the resulting contours:

The input values for this graphic are strings of the sort one normally uses to express complex numbers.

The intersections of the two sets of contours represent the solutions to the defining equation on the complex plane as a function of the two complex input parameters. The contours for the two special cases of the ordinary Lambert function can be found by simply setting each parameter to zero in turn.

In general there will be a line of one or more intersections in the middle of the graphic. Above and below this set of solutions will be pairs of intersections to either side of the middle. As the input parameters are varied, these pairs of intersections move relative to each other, sometimes occurring with the same imaginary part, which will have significant implications for the complex structure of the function. This behavior differs greatly from the ordinary Lambert function, which has only one set of intersections.

Finding the intersections by numerical inversion of the defining equation of course requires a relatively accurate starting point. As a visual experiment, add the first few values of the two special cases to the last graphic. The first special case will be plotted as purple dots and the second special case, with the negative argument, with green dots. In both cases the opacity of the dots decreases with the index value:

Since the defining equation is more complicated that that of the ordinary Lambert function, it is not at all expected that the special cases would give decent approximations to the full equation, but the visual evidence is plain. The reason for including opacity is to show that the values of the index are reversed in sign between the two special cases, a reflection of the difference of sign in the two exponential terms.

When this approximation is good, which of the two special cases applies is determined by the sign of the real part of the Lambert function. The exponential multiplying the parameter x increases in the negative direction and predominates there, while the exponential multiplying the parameter y increases in the positive direction and predominates there.

The approximation is accurate for a variety of inputs, especially on higher branches. Setting the inputs to values with large real or imaginary parts indicates that it is not so accurate in that regime, with the estimates well off the intersections. A second approximation is needed.

When input values are large in absolute value, the right-hand side of the defining equation becomes large enough to make the left-hand side less relevant. One can then write

$\begin{array}{c} x e^{- W (x, y)} + y e^{W (x, y)} \approx 0 \\ e^{2 W (x, y)} \approx \frac{x}{- y} \\ W (x, y) \approx ln \frac{\sqrt{x}}{\sqrt{- y}} = ln \sqrt{x} - ln \sqrt{- y} \end{array}$

Since the ordinary Lambert function can be approximated by a logarithm, this form can be thought of as the sum of the two special cases with an additional factor of two. That motivates attaching the minus sign to y rather than x, which can be verified with an alternate derivation. First let $W = ln Z$ in the defining equation,

$ln Z = \frac{x}{Z} + y Z$

and take partial derivatives of each side:

$\frac{Z_{x}}{Z} = \frac{1}{Z} - \frac{x}{Z^{2}} Z_{x} + y Z_{x} \to Z_{x} = \frac{1}{1 + \frac{x}{Z} - y Z}$

$\frac{Z_{y}}{Z} = - \frac{x}{Z^{2}} Z_{y} + Z + y Z_{y} \to Z_{y} = \frac{Z^{2}}{1 + \frac{x}{Z} - y Z}$

Since both of the denominators are equal, one notices easily that

$Z_{y} = Z^{2} Z_{x}$

which is an exact statement and more importantly separable. Setting $Z (x, y) = f (x) g (y)$ leads to

$f f_{x} = \frac{g_{y}}{g^{3}} = C \to \begin{array}{c} f = \sqrt{2 C x} \\ g = \sqrt{\frac{1}{- 2 C y}} \end{array}$

so that the common constant will cancel out in their product. This leads to

$W (x, y) \approx ln \frac{\sqrt{x}}{\sqrt{- y}} = ln \sqrt{x} - ln \sqrt{- y}$

This approximation is not an exact solution to the defining equation. Presumably it represents the envelope of the exact solutions to the two partial derivatives, a phenomenon that occurs regularly in the context of nonlinear partial differential equations.

Due to the square root, the imaginary part of this approximation changes by $π i$ between branches. With the principal branch indicated by a pink dot and the other branches by gray dots, here is where the values appear relative to the real and imaginary contours, as well as the previous approximation:

In this regime of larger input values this approximation is clearly more accurate than the ordinary Lambert functions as expected. Curiously enough, it holds nicely for the principal branch, apart from cases where either parameter is very small in absolute value, even though it was derived as an asymptotic approximation.

Since this approximation makes the right-hand side of the defining equation identically zero, the generalized Lambert function must contain an additional term for its defining equation to hold. This means that exactly

$W (x, y) = ln \frac{\sqrt{x}}{\sqrt{- y}} + u (x, y)$

where the additional term is expected to be small when the asymptotic approximation is good. Putting this expression in the defining equation gives

$ln \frac{\sqrt{x}}{\sqrt{- y}} + u = - 2 \sqrt{x} \sqrt{- y} sinh u$

which is a hyperbolic form of the Kepler equation. As an extremely interesting aside, that means the Kepler equation is exactly soluble by the generalized Lambert function introduced here.

For small u, keeping only the first-order term of the expansion of the right-hand side of the last equation gives

$u (x, y) \approx \frac{ln \sqrt{x} - ln \sqrt{- y}}{1 + 2 \sqrt{x} \sqrt{- y}}$

on the principal branch. This added term vanishes for large values of either independent variable, since a positive power eventually grows more quickly than a logarithm. Applying l’Hôpital’s rule to the determining factors in this particular case,

$lim_{x \to \infty} \frac{ln x}{\sqrt{x}} = lim_{x \to \infty} \frac{1}{x} \times \frac{1}{2 x^{3 / 2}} = 0$

and likewise for the factors in y. This means that this asymptotic approximation becomes exact asymptotically. This behavior differs considerably from the ordinary Lambert function, which asymptotically requires beyond a simple logarithm an additional term of the logarithm of a logarithm.

For higher branches include the behavior of the logarithm, taking account of the factor of one-half from the square root:

$W_{n} (x, y) = ln \frac{\sqrt{x}}{\sqrt{- y}} + n π i + u (x, y)$

Putting this in the defining equation gives exactly

$ln \frac{\sqrt{x}}{\sqrt{- y}} + u + n π i = - 2 (- 1)^{n} \sqrt{x} \sqrt{- y} sinh u$

and the same first-order expansion leads to the expression

$u_{n} (x, y) \approx \frac{ln \sqrt{x} - ln \sqrt{- y} \pm n π i}{1 + 2 (- 1)^{n} \sqrt{x} \sqrt{- y}}$

which vanishes for large values of the independent variables for the same reasons. This expression can further be used as a guideline for when to switch from the asymptotic approximation to the ordinary Lambert functions: the latter will be more efficient when this value exceeds some small value of about unity.

Knowing that the asymptotic approximation becomes exact asymptotically, it can be used to order the sheets of this complex function, just as one does for the logarithm and ordinary Lambert function. The question remains as to how to choose which of the two Lambert functions forms a continuous Riemann sheet with the asymptotic approximation. A naive approach, that actually turns out to be useful, is to calculate the asymptotic approximation along with both ordinary Lambert functions, then return whichever of these latter two has an imaginary part nearer to that of the asymptotic approximation that labels the sheet.

For this procedure the starting point for numerical inversion has this appearance on the complex plane as a function of x,

where the colors of sections correspond to the colored dots appearing in previous graphics. With this starting point, the individual branches of the generalized Lambert function as a function of x look like this:

The options to include adjacent sheets let one see immediately how the sheets connect to one another. The coloring here is for distinct contrast.

Since the asymptotic approximation is a logarithm, the branch cut on the negative real axis is expected. As one zooms out by increasing the range of the graphic, the function clearly approaches the form of a simple logarithm.

What is not immediately expected is the appearance of a second branch cut that moves dynamically when the input parameters are altered. While its presence can be inferred from the observation above that intersections of contours representing solutions sometimes have the same imaginary part, seeing it directly is somewhat dramatic. This dynamic branch cut connects successive pairs of sheets in both the real and imaginary parts. It is rather unusual in that the real part has a gap across the branch cut while the imaginary parts intersect.

In order to make the dynamic branch cut more evident, a coloring scheme by argument is useful, but not by the argument itself. For higher branches in either direction the argument quickly approaches a limit, making coloring by argument excessively monochromatic: the same phenomenon occurs for both the logarithm and the ordinary Lambert function. Instead the coloring will be relative to the phase of a set point on the real axis, say unity, and scaled to produce the most dramatic effect.

With coloring as described the previous graphic becomes

The nature of the dynamic branch cut is now clearly seen by the change in coloring across it.

For completeness the last two graphics will be repeated as a function of y. First the monochromatic sheets,

and then the relative coloring by argument:

The behavior of the function is essentially the same, apart from the main logarithmic branch cut being located on the positive real axis due to the difference in sign of the independent variable.

Close inspection of the interactive graphics depicting individual sheets of the generalized Lambert function might show gaps or irregularities. These are points where the numerical inversion jumps to incorrect sheets and could stand improvement. For the purpose of an overall visualization of the behavior of the function, however, these will be ignored for the time being.

The purpose of this presentation is to understand the nature of this function well enough to exhibit this overall visualization, and this has been achieved. A deeper question arises in this context: how does one characterize mathematically the dynamic branch cut revealed? This appears to be far from a trivial question and will be left for a future presentation.

And at least one further task remains: a good name for this function. How about the double Lambert function?